knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, 
                      fig.height = 4, fig.width = 7)

library(terra)
library(sf)
library(data.table)
library(tidyverse)
library(here)
library(oharac)
source(here('common_fxns.R'))

Summary

Using functional entities defined in a prior script based on a suite of categorical traits, in combination with spatial distribution of species from AquaMaps and IUCN, calculate functional diversity metrics based on Mouillot et al (2014) and code from Sebastien Villeger: functional richness, functional vulnerability, functional redundancy, and functional over-redundancy.

Methods

Summarize functional entity membership per cell across IUCN and AquaMaps

Included species are:

  • coded for functional entity membership
  • mapped in one (or both) of the species distribution mapsets
  • included the trait-based vulnerability set

Inspect bin membership for FE traits

spp_info <- assemble_spp_info_df() %>%
  select(species, taxon, fe_id, mob, trp, len, wcol) %>%
  distinct() %>%
  mutate(mob = factor(mob, levels = c('ses', 'sed', 'mob', 'mig')),
         trp = factor(trp, levels = c('pri', 'sec', 'ter', 'apex')),
         len = factor(len, levels = c('tiny', 'vsm', 'sm', 'med', 'ml', 'lg', 'vlg', 'huge')))

n_fes_possible <- n_distinct(spp_info$mob) * 
  n_distinct(spp_info$trp) * 
  n_distinct(spp_info$len) * 
  n_distinct(spp_info$wcol)
n_fes_realized <- spp_info$fe_id %>% n_distinct()

Overall, with this set of traits used to define functional entities, we could have up to 512 FEs, only 339 (66.2%) of which are realized in this set of species. Note that the Mouillot 2014 paper includes two additional traits bringing the total number of possible FEs up to 5670, only 646 of which are realized by the 6316 tropical fish spp in their study.

Adult body length

knitr::kable(t(table(spp_info$len)))
tiny vsm sm med ml lg vlg huge
2174 5959 4179 3595 2555 1370 906 406

Trophic level

knitr::kable(t(table(spp_info$trp)))
pri sec ter apex
1647 4247 12810 2440

Adult mobility

knitr::kable(t(table(spp_info$mob)))
ses sed mob mig
2539 4378 10709 3518

Water col position

knitr::kable(t(table(spp_info$wcol)))
ben bp pel rf
12702 1130 2233 5079

Load species ranges and summarize to functional entities per cell

Here we use the IUCN species range maps rasterized to 10 km Mollweide, and AquaMaps species distributions reprojected to 10 km Mollweide. Priority: IUCN species maps, then AquaMaps with occur_cells ≥ 10. These map filenames are already set up in the assemble_spp_info_df() function.

Read in functional entity assignments and join to spatial ranges. Per cell per FE, calculate the number of species for later calculation into functional richness, redundancy, vulnerability. Break the process into chunks to avoid overloading memory.

process_spp_to_fe <- function(fe, spp_maps) { ### fe <- fe_vec[1]
  dt_out <- spp_maps %>%
    filter(fe_id == fe) %>%
    data.table() %>%
    .[ , .(n_spp_fe = length(unique(species))),
       by = .(cell_id, fe_id)]
  return(dt_out)
}
fe_summary_file <- here_anx('func_entities', 'fe_cell_summary.csv')
# unlink(fe_summary_file)
if(!file.exists(fe_summary_file)) {

  spp_info_df <- assemble_spp_info_df(fe_only = TRUE, vuln_only = TRUE)
  ### table(spp_info_df %>% select(species, len) %>% distinct() %>% .$len)
  # huge   lg  med   ml   sm tiny  vlg  vsm 
  #  406 1370 3595 2555 4179 2174  906 5959 

  ### Chunk out and save smaller files to avoid loading in EVERY spp map
  ### all at once!
  n_chunks <- 40
  spp_to_chunks <- spp_info_df %>%
    select(species) %>%
    distinct() %>%
    mutate(chunk = rep(1:n_chunks, length.out = n()))
  
  spp_clean <- spp_info_df %>%
    select(species, fe_id, map_f, taxon) %>%
    distinct() %>%
    left_join(spp_to_chunks, by = 'species')
  
  spp_fes <- spp_clean %>%
    select(species, fe_id) %>%
    distinct()

  tmp_chunk_filestem <- here('tmp/fe_summary_chunk_%s.csv')
  for(i in 1:n_chunks) { ### i <- 1
    tmp_chunk_f <- sprintf(tmp_chunk_filestem, i)
    if(file.exists(tmp_chunk_f)) {
      message('Temp file ', basename(tmp_chunk_f), ' exists... skipping!')
      next()
    }
    message('Assembling species maps and binding FEs for chunk ', i, ' of ', n_chunks, '...')
    spp_chunk <- spp_clean %>% filter(chunk == i)
    chunk_spp_maps <- collect_spp_rangemaps(spp_vec  = spp_chunk$species, 
                                            file_vec = spp_chunk$map_f,
                                            idcol = 'species',
                                            ) %>%
      dt_join(spp_fes, by = 'species', type = 'left') %>%
      distinct()
    
    fe_vec <- chunk_spp_maps$fe_id %>% unique() %>% sort()
    
    message('... summarizing FR map for chunk ', i, ': ', length(fe_vec), ' FEs in ', 
      nrow(chunk_spp_maps), ' cell observations...')
    chunk_fe_sumlist <- parallel::mclapply(fe_vec, mc.cores = 20,
                                           FUN = process_spp_to_fe,
                                           spp_maps = chunk_spp_maps)
    if(check_tryerror(chunk_fe_sumlist)) {
      stop('Try-errors detected in FE summary loop!')
    }
    
    message('... binding chunk and writing chunk summary to ', basename(tmp_chunk_f), '...')
    chunk_fe_summary <- chunk_fe_sumlist %>% data.table::rbindlist()
    fwrite(chunk_fe_summary, tmp_chunk_f)
  }
  
  message('... reading temp chunks, binding, and summarizing by fe and cell...')
  fe_tmp_fs <- list.files(dirname(tmp_chunk_filestem), 
                          pattern = 'fe_summary_chunk', full.names = TRUE)
  fe_summary_list <- parallel::mclapply(fe_tmp_fs, FUN = fread, mc.cores = 10)
  fe_summary <- rbindlist(fe_summary_list) %>%
    .[ , .(n_spp_fe = sum(n_spp_fe)),
       by = .(cell_id, fe_id)]
  
  message('... writing summary to ', basename(fe_summary_file), '...')
  fwrite(fe_summary, fe_summary_file)
  ### unlink(fe_tmp_fs)
}

Calculate functional diversity metrics

Calculate functional diversity metrics from IUCN and AquaMaps summaries of functional entity membership per cell.

This set of calculations replicates the results of Villeger’s function but without resorting to a presence/absence matrix, instead just relying on the original AquaMaps species-cell table.

  • spp_fe_df joined to am_spp_cell;
  • group by cell and functional group;
  • calculate number of spp per functional group;
  • group by cell;
  • calculate functional entity metrics -
    • Species richness
    • Functional entity richness
    • Functional vulnerability (% of FEs represented by a single spp)
    • Weighted functional vulnerability (new: vuln of an FE is \(\frac{1}{2^{n-1}}\), take average over all FE)
    • Functional redundancy (mean spp count per FE)
    • Functional over-redundancy (sum of spp per FE above mean, divided by total number of spp)
fe_metrics_file <- here_anx('func_entities/fe_metrics_per_cell.csv')

### unlink(fe_metrics_file)
if(!file.exists(fe_metrics_file)) {
  # system.time({
  fe_sum_total <- data.table::fread(fe_summary_file)

  ### there are 6.5 million cells.  Chop into 500k instances and mclapply to summarize
  # cell_id_vec <- 1:ncell(raster::raster(here('_spatial/ocean_area_mol.tif')))
  chunk_size <- 250000
  n_chunks <- 6.5e6 / chunk_size

  nspp_tmp_f <- here('tmp/nspp_df_tmp.csv')
  if(!file.exists(nspp_tmp_f)) {
    message('Calculating species richness map...')
    nspp_list <- parallel::mclapply(1:n_chunks, mc.cores = ceiling(n_chunks / 4),
                   FUN = function(n) { ### n <- 1
                     message('... spp richness chunk ', n, ' of ', n_chunks, '...')
                     cell_id_min <- (n - 1) * chunk_size + 1
                     cell_id_max <- n * chunk_size
                     nspp_sum <- fe_sum_total %>%
                       filter(between(cell_id, cell_id_min, cell_id_max)) %>%
                       data.table() %>%
                       .[ , .(n_spp = sum(n_spp_fe, na.rm = TRUE),
                              n_fe  = n_distinct(fe_id)),
                          by = 'cell_id']
                     }) 
    
    if(check_tryerror(nspp_list)) stop('Try errors detected!')
    nspp_df <- nspp_list %>%
      data.table::rbindlist()
    
    write_csv(nspp_df, nspp_tmp_f)
  } else { 
    message('File exists: ', basename(nspp_tmp_f))
    nspp_df <- data.table::fread(nspp_tmp_f)
  }
  
  fvuln_tmp_f <- here('tmp/fvuln_df_tmp.csv')
  if(!file.exists(fvuln_tmp_f)) {

    message('Calculating functional vulnerability map...')
    fv_list <- parallel::mclapply(1:n_chunks, mc.cores = ceiling(n_chunks / 4),
                   FUN = function(n) { ### n <- 1
                     message('... funct vuln chunk ', n, ' of ', n_chunks, '...')
                     cell_id_min <- (n - 1) * chunk_size + 1
                     cell_id_max <- n * chunk_size
                     f_vuln <- fe_sum_total %>%
                       filter(between(cell_id, cell_id_min, cell_id_max)) %>%
                       oharac::dt_join(nspp_df, by = 'cell_id', type = 'left') %>%
                       data.table() %>%
                       .[ , .(f_vuln  = sum(n_spp_fe == 1) / n_fe, 
                              f_wvuln = mean(0.5^(n_spp_fe - 1)),
                              f_red   = mean(n_spp_fe)),
                       by = .(cell_id, n_fe, n_spp)]
                   }) 
    if(check_tryerror(fv_list)) stop('Try errors detected!')
    f_vuln_df <- fv_list %>%
      data.table::rbindlist()
    write_csv(f_vuln_df, fvuln_tmp_f)
  } else { 
    message('File exists: ', basename(fvuln_tmp_f))
    f_vuln_df <- data.table::fread(fvuln_tmp_f)
  }
  

  fred_tmp_f <- here('tmp/f_red_df_tmp.csv')
  if(!file.exists(fred_tmp_f)) {
    message('Calculating functional overredundancy map...')
    f_overred_list <- parallel::mclapply(1:n_chunks, mc.cores = ceiling(n_chunks / 4),
                FUN = function(n) { ### n <- 1
                  message('... funct overredundancy chunk ', n, ' of ', n_chunks, '...')
                  cell_id_min <- (n - 1) * chunk_size + 1
                  cell_id_max <- n * chunk_size
                  f_overred <- fe_sum_total %>%
                    filter(between(cell_id, cell_id_min, cell_id_max)) %>%
                    oharac::dt_join(f_vuln_df, by = 'cell_id', type = 'left') %>%
                    data.table() %>%
                    .[ , f_over := ifelse(n_spp_fe > f_red, n_spp_fe - f_red, 0)] %>%
                    .[ , .(f_overred = sum(f_over) / sum(n_spp_fe)),
                       by = 'cell_id']
                })
    if(check_tryerror(f_overred_list)) stop('Try errors detected!')
    f_overred_df <- f_overred_list %>%
      data.table::rbindlist()
    write_csv(f_overred_df, fred_tmp_f)
  } else { 
    message('File exists: ', basename(fred_tmp_f))
    f_overred_df <- data.table::fread(fred_tmp_f)
  }

  message('Combining functional diversity metrics maps...')
  fe_metrics_df <- f_vuln_df %>%
    oharac::dt_join(f_overred_df, by = 'cell_id', type = 'full')
  
  message('Combined set: ', nrow(fe_metrics_df), ' rows!')
  write_csv(fe_metrics_df, fe_metrics_file)
}

Map functional diversity metrics

Create and save rasters of functional diversity metrics

fe_metrics_df <- data.table::fread(fe_metrics_file)
### note that while there are 3,684,273 cells in this dataframe,
### 4273 of those are the Caspian Sea, so when those are dropped,
### it results in 3,680,000 exactly... ugh, weird

n_spp_rast   <- map_to_mol(fe_metrics_df, which = 'n_spp')
n_fe_rast    <- map_to_mol(fe_metrics_df, which = 'n_fe')
f_vuln_rast  <- map_to_mol(fe_metrics_df, which = 'f_vuln')
f_wvuln_rast <- map_to_mol(fe_metrics_df, which = 'f_wvuln')
f_red_rast   <- map_to_mol(fe_metrics_df, which = 'f_red')
f_overred_rast <- map_to_mol(fe_metrics_df, which = 'f_overred')

writeRaster(n_spp_rast,   here('_output/func_entities/n_spp.tif'),   
            overwrite = TRUE)
writeRaster(n_fe_rast,    here('_output/func_entities/n_fe.tif'),    
            overwrite = TRUE)
writeRaster(f_vuln_rast,  here('_output/func_entities/f_vuln.tif'), 
            overwrite = TRUE)
writeRaster(f_wvuln_rast, here('_output/func_entities/f_wvuln.tif'), 
            overwrite = TRUE)
writeRaster(f_red_rast,   here('_output/func_entities/f_red.tif'),   
            overwrite = TRUE)
writeRaster(f_overred_rast, here('_output/func_entities/f_overred.tif'), 
            overwrite = TRUE)

Examine distributions of each raster.

Plot each raster. NOTE: These no longer include spp that are missing from vulnerability calculations - so the maps in _output/nspp_maps (vuln \(\cap\) FE) should be identical to these in output/func_entities.